
rm(list=ls())
gc()
setwd("/hpc/group/pbenfeylab/CheWei/")


as.numeric(system("awk '/MemFree/ {print $2}' /proc/meminfo", intern=TRUE))


suppressMessages(library(Seurat))
suppressMessages(library(tricycle))
suppressMessages(library(cowplot))
suppressMessages(library(ggplot2))
suppressMessages(library(scattermore))
suppressMessages(library(scater))
suppressMessages(library(cowplot))
suppressMessages(library(RColorBrewer))
suppressMessages(library(grid))
suppressMessages(library(gplots))
suppressMessages(library(tidyverse))

prep <- function(seu){
seu$time.anno.Li.crude <- gsub("_.*$","",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Columella","Distal Columella_Columella",seu$time.celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Distal Lateral Root Cap","Distal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Columella","Proximal Columella_Columella",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("Proximal Lateral Root Cap","Proximal Lateral Root Cap_Lateral Root Cap",seu$celltype.anno.Li.crude)
seu$celltype.anno.Li.crude <- gsub("^.*_","",seu$celltype.anno.Li.crude)
    return(seu)
}
plot_anno <- function(rc.integrated){
order <- c("Quiescent Center", "Ground Tissue","Columella", "Lateral Root Cap", "Atrichoblast", "Trichoblast", "Cortex", "Endodermis", "Phloem","Protophloem", "Xylem", "Procambium","Pericycle","Phloem Pole Pericycle", "Protoxylem", "Metaxylem", "Unknown")
palette <- c("#9400d3", "#DCD0FF","#5ab953", "#bfef45", "#008080", "#21B6A8", "#82b6ff", "#0000FF","#e6194b", "#dd77ec", "#9a6324", "#ffe119", "#ff9900", "#ffd4e3", "#9a6324", "#ddaa6f", "#EEEEEE")
rc.integrated$celltype.anno.Li.crude <- factor(rc.integrated$celltype.anno.Li.crude, levels = order[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]) 
color <- palette[sort(match(unique(rc.integrated$celltype.anno.Li.crude),order))]
p1 <- DimPlot(rc.integrated, reduction = "umap", group.by = "celltype.anno.Li.crude", cols=color)
p2 <- DimPlot(rc.integrated, reduction = "umap", group.by = "time.anno.Li.crude", order = c("Distal Columella","Proximal Columella","Distal Lateral Root Cap","Proximal Lateral Root Cap","Maturation","Elongation", "Transition Domain", "Proliferation Domain"),
        cols = c('#F7E7B0','#FFC400','#2B871F','#005E3B',"#deebf7", "#3182bd", '#fee0d2','#de2d26'))
p3 <- DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.time.group", cols=brewer.pal(10,"Spectral"))
p4 <- DimPlot(rc.integrated, reduction = "umap", group.by = "branch.anno")
options(repr.plot.width=16, repr.plot.height=12)
gl <- lapply(list(p1, p2, p3, p4), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}


sessionInfo()

R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: AlmaLinux 9.3 (Shamrock Pampas Cat)

Matrix products: default
BLAS/LAPACK: /hpc/group/pbenfeylab/ch416/miniconda3/envs/seu4/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] forcats_0.5.2               stringr_1.5.1              
 [3] dplyr_1.1.3                 purrr_1.0.2                
 [5] readr_2.1.3                 tidyr_1.3.0                
 [7] tibble_3.2.1                tidyverse_1.3.2            
 [9] gplots_3.1.3                RColorBrewer_1.1-3         
[11] scater_1.26.1               scuttle_1.8.0              
[13] scattermore_1.2             ggplot2_3.4.4              
[15] cowplot_1.1.1               tricycle_1.6.0             
[17] SingleCellExperiment_1.20.0 SummarizedExperiment_1.28.0
[19] Biobase_2.58.0              GenomicRanges_1.50.0       
[21] GenomeInfoDb_1.34.8         IRanges_2.32.0             
[23] S4Vectors_0.36.0            BiocGenerics_0.44.0        
[25] MatrixGenerics_1.10.0       matrixStats_1.1.0          
[27] SeuratObject_4.1.3          Seurat_4.1.1.9001          

loaded via a namespace (and not attached):
  [1] utf8_1.2.4                spatstat.explore_3.2-5   
  [3] reticulate_1.34.0         tidyselect_1.2.0         
  [5] RSQLite_2.3.1             AnnotationDbi_1.60.0     
  [7] htmlwidgets_1.6.2         BiocParallel_1.32.5      
  [9] Rtsne_0.16                munsell_0.5.0            
 [11] ScaledMatrix_1.6.0        codetools_0.2-19         
 [13] ica_1.0-3                 pbdZMQ_0.3-8             
 [15] future_1.33.0             miniUI_0.1.1.1           
 [17] withr_2.5.2               spatstat.random_3.2-1    
 [19] colorspace_2.1-0          progressr_0.14.0         
 [21] uuid_1.1-0                ROCR_1.0-11              
 [23] tensor_1.5                listenv_0.9.0            
 [25] repr_1.1.4                GenomeInfoDbData_1.2.9   
 [27] polyclip_1.10-6           bit64_4.0.5              
 [29] parallelly_1.36.0         vctrs_0.6.4              
 [31] generics_0.1.3            circular_0.4-95          
 [33] timechange_0.1.1          R6_2.5.1                 
 [35] ggbeeswarm_0.7.1          rsvd_1.0.5               
 [37] bitops_1.0-7              spatstat.utils_3.0-4     
 [39] cachem_1.0.8              DelayedArray_0.24.0      
 [41] assertthat_0.2.1          promises_1.2.1           
 [43] scales_1.2.1              googlesheets4_1.0.1      
 [45] beeswarm_0.4.0            gtable_0.3.4             
 [47] beachmat_2.14.0           globals_0.16.2           
 [49] goftest_1.2-3             rlang_1.1.2              
 [51] splines_4.2.2             lazyeval_0.2.2           
 [53] gargle_1.2.1              broom_1.0.2              
 [55] spatstat.geom_3.2-7       modelr_0.1.10            
 [57] reshape2_1.4.4            abind_1.4-5              
 [59] backports_1.4.1           httpuv_1.6.12            
 [61] tools_4.2.2               ellipsis_0.3.2           
 [63] ggridges_0.5.4            Rcpp_1.0.11              
 [65] plyr_1.8.9                base64enc_0.1-3          
 [67] sparseMatrixStats_1.10.0  zlibbioc_1.44.0          
 [69] RCurl_1.98-1.6            deldir_1.0-9             
 [71] pbapply_1.7-2             viridis_0.6.4            
 [73] zoo_1.8-12                haven_2.5.1              
 [75] ggrepel_0.9.4             cluster_2.1.4            
 [77] fs_1.6.3                  magrittr_2.0.3           
 [79] data.table_1.14.8         RSpectra_0.16-1          
 [81] reprex_2.0.2              lmtest_0.9-40            
 [83] RANN_2.6.1                googledrive_2.0.0        
 [85] mvtnorm_1.1-3             ggnewscale_0.4.8         
 [87] fitdistrplus_1.1-11       hms_1.1.2                
 [89] patchwork_1.1.3           mime_0.12                
 [91] evaluate_0.23             xtable_1.8-4             
 [93] readxl_1.4.1              fastDummies_1.7.3        
 [95] gridExtra_2.3             compiler_4.2.2           
 [97] KernSmooth_2.23-20        crayon_1.5.2             
 [99] htmltools_0.5.7           tzdb_0.3.0               
[101] later_1.3.1               lubridate_1.9.0          
[103] DBI_1.1.3                 dbplyr_2.2.1             
[105] MASS_7.3-58.3             boot_1.3-28.1            
[107] Matrix_1.5-4              cli_3.6.1                
[109] parallel_4.2.2            igraph_1.5.1             
[111] pkgconfig_2.0.3           sp_2.1-1                 
[113] IRdisplay_1.1             plotly_4.10.3            
[115] spatstat.sparse_3.0-3     xml2_1.3.3               
[117] vipor_0.4.5               XVector_0.38.0           
[119] rvest_1.0.3               digest_0.6.33            
[121] sctransform_0.4.1         RcppAnnoy_0.0.21         
[123] spatstat.data_3.0-3       Biostrings_2.66.0        
[125] cellranger_1.1.0          leiden_0.4.3             
[127] uwot_0.1.16               DelayedMatrixStats_1.20.0
[129] shiny_1.7.5.1             gtools_3.9.4             
[131] lifecycle_1.0.4           nlme_3.1-162             
[133] jsonlite_1.8.7            BiocNeighbors_1.16.0     
[135] viridisLite_0.4.2         fansi_1.0.5              
[137] pillar_1.9.0              lattice_0.21-8           
[139] KEGGREST_1.38.0           fastmap_1.1.1            
[141] httr_1.4.7                survival_3.4-0           
[143] glue_1.6.2                png_0.1-8                
[145] bit_4.0.5                 stringi_1.8.1            
[147] blob_1.2.4                RcppHNSW_0.5.0           
[149] BiocSingular_1.14.0       caTools_1.18.2           
[151] memoise_2.0.1             IRkernel_1.3.1.9000      
[153] irlba_2.3.5.1             future.apply_1.11.0


new.ref <- readRDS("./tradeseq/Tricycle_Reference_Arabidopsis_Root.rds")


run_tricycle <- function(sample){
    seu <- readRDS(paste0("./COPILOT_other_projects/to_be_integrate/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}


sample.list <- c("sc_283","sc_284","sc_285","sc_286","sc_287","sc_288","sc_289","sc_290","sc_291","sc_292","sc_293","sc_294")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


run_tricycle <- function(sample){
    seu <- readRDS(paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT.rds"))
    sce <- as.SingleCellExperiment(seu)
    
    ## Tricycle estimation
    sce <- project_cycle_space(sce,ref.m = new.ref)
    sce <- estimate_cycle_position(sce)
    
    sce$tricyclePosition <- (sce$tricyclePosition/pi) + 1.8
    sce$tricyclePosition[which(sce$tricyclePosition>=2)] <- sce$tricyclePosition[which(sce$tricyclePosition>=2)] - 2
    sce$tricyclePosition <- sce$tricyclePosition*pi
    
    ## Plot
    print(plot_ccposition_den(sce$tricyclePosition,
                    sce$orig.ident, 'sample',
                    bw = 10, fig.title = "Kernel density of \u03b8") +
      theme_bw(base_size = 14))
    
    ## Assign position to stage
    tricycleCCStage <- sce$tricyclePosition/pi

    G1_idx <- which(tricycleCCStage < 1 & tricycleCCStage >= 0)
    S_idx <- which(tricycleCCStage >= 1 & tricycleCCStage < 1.55)
    G2_idx <- which(tricycleCCStage >= 1.55 & tricycleCCStage <= 2)
                    
    tricycleCCStage[G1_idx]="G1/G0"
    tricycleCCStage[S_idx]="S"
    tricycleCCStage[G2_idx]="G2M"
    
    sce$tricycleCCStage <- tricycleCCStage
    
    ## Schwabe Stage
    #sce <- estimate_Schwabe_stage(sce,cycleGene.l = cc.gene.list)
    
    ## Table
    #seu$time.anno.Li <- gsub("_.*$","",seu$time.celltype.anno.Li)
    #print(table(seu$time.anno.Li, sce$tricycleCCStage))
    
    ## Prepare seu
    #seu <- prep(seu)
    #plot_anno(seu)
    
    ## Sce to seu
    seu$tricyclePosition <- sce$tricyclePosition/pi
    seu$tricycleCCStage <- sce$tricycleCCStage
    #seu$SchwabeCCStage <- sce$CCStage
    
    ## Tricycle position on umap
    p <- plot_emb_circle_scale(sce, dimred = "UMAP",
                           point.size = 1, point.alpha = 0.9) +
      theme_void(base_size = 14)
    legend <- circle_scale_legend(text.size = 5, alpha = 0.9)
    print(plot_grid(p, legend, ncol = 2, rel_widths = c(1, 0.5)))
    
    ## Write annotations as csv
    write.csv(seu@meta.data[,c("tricyclePosition", "tricycleCCStage")], paste0("./tricycle/",sample,".csv"), quote=FALSE)
    #saveRDS(seu,paste0("./scRNA-seq/Seurat_Objects/",sample,"_COPILOT_seu4.rds"))
}


sample.list <- c("sc_155","sc_156","sc_157","sc_158","sc_159","sc_160","sc_161","sc_162","sc_163","sc_164","sc_165","sc_166","sc_167","sc_168","sc_169","sc_170","sc_171","sc_172","sc_173","sc_174","sc_175","sc_176","sc_177","sc_178","sc_179","sc_180")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 240.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”


sample.list <- c("sc_57","sc_58","sc_59","sc_60","sc_61","sc_62","sc_63","sc_64","sc_66")


for (i in 1:length(sample.list)){
 run_tricycle(sample.list[i])
}

The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”
The number of projection genes found in the new data is 245.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 243.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 242.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 244.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

The number of projection genes found in the new data is 223.

Warning message in brewer.pal(nlevels(color_var.v), "Set1"):
“minimal value for n is 3, returning requested palette with 3 different levels
”

	used	(Mb)	gc trigger	(Mb)	max used	(Mb)
Ncells	625269	33.4	1361470	72.8	1142193	61.0
Vcells	1159062	8.9	8388608	64.0	1802279	13.8

Load reference¶

For sample under ./scRNA-seq/Seurat_Objects/¶